
#include <stdio.h>
#include <math.h>
#include <time.h>

#define pi 3.1415926535877932384626
#define N 3500
//#define G 6.673e-11
#define G 0.001
//#define theta i*2*pi/N

/** initial mass, velocity and position */
void
initial (double m[], double v[][2], double p[][2])
{
  int i;
  for (i = 0; i < N; i++)
    {
      m[i] = 1;
      double theta = i * 2 * pi / N;
      // distributed around unit circle
      p[i][0] = cos (theta);
      p[i][1] = sin (theta);
      // velocity point to tangent direction
      //v[i][0] = -sin(theta);
      //v[i][1] = cos(theta);
      // velocity point to origin, or normal direction
      v[i][0] = -p[i][0];
      v[i][1] = -p[i][1];

    }
}

int
main ()
{
  static double m[N];
  static double v[N][2];
  static double p[N][2];
  static double f[N][N][2];
  static double ftot[N][2];
  static double a[N][2];

  int i, j, k;
  double dt = 0.0001;
  double mv0 = 0, mv1 = 0, mvr0 = 0, mvr1 = 0;

  initial (m, v, p);
  /*
     for (i=0;i<N;i++){
     double r = sqrt(p[i][0]*p[i][0]+p[i][1]*p[i][1]);
     mv0 +=v[i][0];
     mv1 +=v[i][1];
     mvr0 +=v[i][0]*r;
     mvr1 += v[i][1]*r;
     }

     printf("momentum: %f  %f \n angle momentum: %f %f\n", mv0, mv1, mvr0, mvr1); */
  double iter = 40;
  double r_ij0, r_ij1, r_ij, rcube;
  clock_t begin, end;
  double timecost;

  begin = clock ();

  for (k = 0; k < iter; k++)
    {
      //compute f_ij for all elements
      for (i = 0; i < N; i++)
	{
	  for (j = i + 1; j < N; j++)
	    {

	      r_ij0 = p[j][0] - p[i][0];
	      r_ij1 = p[j][1] - p[i][1];
	      r_ij = sqrt (r_ij0 * r_ij0 + r_ij1 * r_ij1);
	      rcube = r_ij * r_ij * r_ij;
	      f[i][j][0] = G * m[i] * m[j] * r_ij0 / rcube;
	      f[i][j][1] = G * m[i] * m[j] * r_ij1 / rcube;

	      f[j][i][0] = -f[i][j][0];
	      f[j][i][1] = -f[i][j][1];
	    }
	}

      //update velocity and position
      for (i = 0; i < N; i++)
	{
	  ftot[i][0] = 0;
	  ftot[i][1] = 0;
	  for (j = 0; j < N; j++)
	    {
	      ftot[i][0] += f[i][j][0];
	      ftot[i][1] += f[i][j][1];
	    }
	  a[i][0] = ftot[i][0] / m[i];
	  a[i][1] = ftot[i][1] / m[i];

	  p[i][0] += v[i][0] * dt;
	  p[i][1] += v[i][1] * dt;

	  v[i][0] += a[i][0] * dt;
	  v[i][1] += a[i][1] * dt;


	}
    }
  /*for (  i=0; i<N; i++){
     //printf("r: %f  position: %f %f\n velocity: %f %f\n acc: %f %f\n", sqrt(p[i][0]*p[i][0]+p[i][1]*p[i][1]), p[i][0], p[i][1], v[i][0],v[i][1], a[i][0], a[i][1]);
     double r = sqrt(p[i][0]*p[i][0]+p[i][1]*p[i][1]);
     mv0 +=v[i][0];
     mv1 +=v[i][1];
     mvr0 +=v[i][0]*r;
     mvr1 += v[i][1]*r;
     }

     printf("momentum: %f  %f \n angle momentum: %f %f\n", mv0, mv1, mvr0, mvr1); */
  end = clock ();
  timecost = (double) (end - begin) / CLOCKS_PER_SEC;
  double t_iter = timecost / iter;
  FILE *fp;
  fp = fopen ("result", "a");
  fprintf (fp, "%d  %f\n", N, t_iter);
  printf ("%d body, t_iter=%f s\n", N, t_iter);
  fclose (fp);
  return 1;

}
